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ABSTRACT 

The origin of the supermassive black holes that power the most distant quasars observed is largely 
unknown. One hypothesis is that they grew rapidly from intermediate-mass seeds (~ 100 M ) left 
by the first stars. However, some previous studies argued that accretion onto these black holes was 
too low to build up the mass due to strong suppression by radiative feedback. Here, we re-exam the 
accretion process of such a black hole embedded in a primordial gas cloud, by considering a wide 
range of physical and numerical parameters not explored before. We find that, while radiative heating 
and pressure indeed suppress accretion effectively, self-gravity of the gas eventually overcomes the 
feedback effects and boosts the accretion to the Eddington rate after one free-fall timescale of the 
cloud. Moreover, for a given black hole mass, there exists a critical density above which the accretion 
can reach Eddington limit. Furthermore, we find a universal correlation between black hole accretion 
rate and ambient gas density, which may serve as a realistic recipe for black hole growth in simulations. 
Subject headings: accretion - black hole physics - radiative feedback - hydrodynamics - methods: 
numerical - quasars: high redshift 



1. INTRODUCTION 

A major recent development in observational cosmol- 
ogy has been the discovery of dozens luminous quasars 
at high redshifts (z > 6) when the universe was less 



than 7% of its current age dFan et al.|200H |2003 Willott 
eT^I][2507l |2010b[ |Mortlock~eTal.||2(Jll| . These quasars 



are believed to be powered by supermassive black holes 
(SMBHs) of - 1O 8_9 M , and the BHs appear to accrete 



_ 2Rl_ 

at near Ed dington rate (e.g., |Fan et al.| 2006| |Willott 



et al.||2010a ). Furthermore, inten se star formation ( |Car 



et al. 



ill! et aygU04{ 

dust (| Walter et 



Walt eretal.|2009|), abundant CO g . 
t al. 1120031 1 Jiang et al.||2006 | IMof 



ms and 



Wang 



2010 ), and super-solar metallicity (iMaiohno et al. 



2005 ) were detected in the quasar hosts, indicating a co 
eval formation of the SMBHs and host galaxies. 

The seeds and growth of these SMBHs, however, are 
unsolved puzzles. It has long been proposed that they 
might grow from remnant s of the very first stars, the so 

|Haim an| 
aTTPOi 



ian| 
U7f 



called PopIII s tars (e.g., jMadau fc Ree s| 
fc Loebl [20011 JVolonteri fc Rees||^005t . 
Volonteri|2010 ) ) . Sophisticated cosmoiogical simulations 
over the past decade suggested that PopIII stars formed 
at early cosmic times z ~ 20 — 50 in minihalos with 
mass — 1O 6 M , an d were likely massive, with mass 
300 M 



30 



2002 

son||2004H ran fc'McK'ee||2 , (j(J41|Cao et a 



V) (e.g. 



Abel et al. 



et ai.l 120081 IBromm et all 120091) ( 



see, 



Bromm & Lar 
2007( 1 Vosmda 



rowever, recent 



ork by |Cfark et al.pOilffOreii' et al][2TiTTl for sugges- 
tions of a smaller mass range) . Massive PopIII stars of 
100— 140 M or above 260 M were predicted to r apidly 
collapse to BHs of - 100 M Q flHeger et aJ.T]|2003| > . The 
subsequent growth of these stellar-mass BHs might play 
an important role in the formation of luminous quasars at 
later times. If a 100 M BH starts to accrete at the Ed- 
dington rate at z > 20, it will attain a final mass greater 
than 10 9 M by z ~ 6, implying a viable explanation for 
the observations of bright quasars in this epoch. 

yuexing@astro.psu.edu 



Black holes are thought to accrete through a disk 
shape d by the outward transfer of the angular momen- 
tum ( Shakura fc Sunyaev| |1973 Rees |1984 Antonucci 
1993 Krolik 1999p . However, m metal-tree, primordial 
galaxies with virial temperature above 10 4 K, conditions 
may exist for the formation of a thick disk of gas at tem- 
perature of 5000 - 10000 K due to insufficie nt cooling 
( |Oh fc Haiman][2002l JVolonteri fe Rees||2005j ). As a re- 
sult, accretion proceeds in a quasi-radial manner onto the 
central BH. Such spherical Bondi-H oyle accretion model 
1944| |Bondi||1952[ ) has been widely em- 



(Bondi & Hoylc 



ployed in cosmo 



ogical simulations on BH growth in the 



early universe (|Li et al.|2007]| Johnson fe Bromm|2007 



Matteo et al.|2008||AlvareTet al. 2009; Sijacki et alT^DO' 



Di Matteo et al 



2011) 



IP] 



In particular, it was shown by 
Li et al.|(|2007p that 100 M BH seeds from PopIII stars 



m gas-rich protogalaxies residing in highly overdense re- 
gions could grow to 10 9 M within 800 million years and 
produce luminous quasars at z — 6 as those observed 



(Li e t al. 2008 Robertson et al. 2007 Narayanan et al. 
2008| ), provided that they accreted at Eddington rates 
for much of their early time. 

There are two primary feedback effects of the radia- 
tion from th e BH that greatly influence the accretion 



Ciotti fc Ostriker|2001||Sazonov et al.|2005 

Alvarez 



process (e.g. 

Ciotti & Ostr'iker|2^07l I Johnson & Bromm|2007 



et al ||2009| |Milosavljevic et al.||2009a|b[ |Park fc Ricott 
2011|)' The first one is the photo-ionization heating of 



the gas within the sonic radius (at which the gas mo- 
tion becomes supersonic), which significantly increases 
the sound speed and reduces the accretion rate. The sec- 
ond one is the radiation pressure from both electron scat- 
tering and photo-ionization, which may drive an outflow 
and halt the accretion completely. To study these feed- 
back processes and their effects on the accretion rate of 
the central BH, one must resolve the spatial scale smaller 
than the sonic radius, which is of the order of 10 14 cm 
(~ 10~ 4 pc) for a 100 M BH. This clearly poses a signif- 
icant challenge for large-scale cosmoiogical simulations. 
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Recent analytical (Milosavljevic et a" 



all and ide- 



alized simu lations (Milosavlj evic et al.|2009b||Park fc Ri 
cotti|2011 ) of accretion onto a 1UU M BH in a gas sphere 



argued that radiation feedback strongly suppressed the 
accretion rate to be at most a small fraction of the Ed- 
dington limit. However, these studies have focused only 
on the regime where the ionization radius is much larger 
than the sonic radius. Furthermore, self-gravity of the 
gas was ignored on the assumption that the BH's grav- 
ity dominates in the ionized zone. But even if the initial 
density is relatively low, as long as the radius of the gas 
sphere is much larger than the ionization region, self- 
gravity would build up the density just outside the ion- 
ization radius to a high level after one free-fall timescale. 

Here, we report new findings of accretion onto a stellar- 
mass BH (A/bh = 100 M Q ) embedded in a primordial gas 
sphere from numerical simulations, which cover a wider 
range of numerical and physical parameters than previ- 
ous studies. Our model includes not only the relevant 
feedback processes, but also self-gravity of the gas. We 
use a modified versio n of the public, grid-bas ed hydro- 



dynamics code VH-1 (Blondin & Lufkin 
the accretion in a spherical geome try. V 



1993} to follow 
TTs based on 



the p iece-wise parabolic method (Colella & Woodward 
1984 1, and uses a Lagrange-Remap approach to solve the 
Euler equations. 

The modifications we made to VH-f include two ma- 
jor aspects: a logarithmic grid to ensure sufficiently high 
resolution near the central hole while covering a large re- 
gion with a reasonable number of grid cells, and a unique 
treatment of radiative feedback by using a finely pre- 
computed two-parameter grid, the ionization parameter 
and gas temperature, to tabulate the cooling, heating, 
and radiation pressure. Our feedbac k algorithm is an 
improvement to those previously used (|Ciotti fe Ostriker 
2001||Sazonov et al.|2005}|Ciotti fe Ostriker|2007| ) in that 



it calculates the cooling, heating, and radiation force 
more accurately. Moreover, it can handle high-density 
clouds accurately and efficiently. These modifications al- 
low us to achieve an unprecedentcdly high spatial reso- 
lution of 10 11 cm f~ 1CP 7 pc), and simulate clouds of 



a wide range of gas density in 10 — 10 cm . In par- 
ticular, the ability to reach densities several orders of 
magnitude higher than the previous limit of 10 7 cm~ 3 is 
critical to our new results. 

The paper is organized as follows. In § 2, we describe 
our analytical considerations and numerical methods in 
detail. In § 3, we first demonstrate that including self- 
gravity enhances the gas density just outside the ioniza- 
tion radius by several orders of magnitude, which leads 
to enhanced accretion rate. We then present various sim- 
ulations with higher but fixed ambient densities without 
including self-gravity, and show that beyond a certain 
critical density, the accretion rate indeed reaches Edding- 
ton limit, and that there exists a correlation between BH 
accretion rate and ambient gas density. We discuss the 
implications of our results to large scale numerical simu- 
lations of SMBH formation, and summarize in § 4. 



2.1. 



2. METHODOLOGY 
Analytical Considerations 



In the analytical work of 
and numerical simulations o 


Milosavl 


evic et al.| (|2009a) 


Milosav 


jevic et ai. 


(2009b) 



and Park & Ricotti (2011), there is an implicit assump- 
tion tnluTtlie^mTzation radius, r; on , is much larger than 
the sonic radius, r s . As we will show later, the conclusion 
that the accretion rate is significantly suppressed by ra- 
diative feedback processes is in fact critically dependent 
on this assumption, more accurately, the assumption that 
''ion S> ^b, where tb is the standard Bondi accretion ra- 
dius ignoring radiative feedback processes: 



rB 



5 x 10 15 — cm. 
Tqa 



Mo 



(1) 



where Mi is the BH mass in unit of 100 M@, and Tq^ is 
the ambient gas temperature at infinity in u nit of 10 4 K. 

the ioniza- 
tion radius r\, 



According to Milosavljevic et al.| (|2009a 
- lon takes the form: 



= 4.7 x 10 



IS 



HII.4.7 



/2/3 & j Orp 
turb 71 * Jl 



2/3, 



,2/3 
HI.3.7 



cm, 



(2) 



where fi on is the average fraction of energy of an ab- 
sorbed photon that goes to photo-ionization, and is 1/3 
for a power-law spectrum with an index of 1.5, as we will 
assume in the present work; I is the BH luminosity rela- 
tive to the Eddington luminosity for Thomson scattering; 
and 715 is the ambient gas density in unit of 10 5 cm -3 . 

By taking other numerical factors to be of order unity, 
we have r ion < r B , if n 5 M 2 > 21 1 / 2 x 10 4 . For a BH 
with Mi = 1, and at near Eddington limit, this means 
n > 2 x 10 9 cm -3 . This density is much higher than 
those co nsid ered in the simulations of Milosav ljevic et al~| 
( |2009b| and |Park fc Ricotti| ( |2011[ ). It is therefore not 
surprising that they found strongly suppressed accre- 
tion rates. However, one can naively imagine that there 
should be qualitative differences in the properties of the 
accretion flow when the r- lon — tb boundary is crossed. 
The reason is that r lon primarily determines the region 
where photon-heating and photo-ionization are impor- 
tant processes, while tb determines the region where the 
BH's gravitational field starts to significantly modify the 
gas flow. When r ion < tb, the radiation feedback is not 
important at radius larger than tb except for Thomson 
scattering. The gas flow can therefore achieve significant 
inward velocity of the order of local sound speed at tb- 
On the other hand, at locations immediately inside n on , 
radiation feedback may slow down the inward velocity 
significantly, creating a strong shocked region between 
rjon and r^. It is this region that may help to maintain 
the accretion rate of the flow at near Eddington limit. 

One critical question is therefore what kind of gas den- 
sity can one expect to find near the accreting BHs. In 
this paper, we suggest that densities as high as, or higher 
than 10 9 cm -3 might be quite common. One reason is 
that the near isothermal collapse of the gas outside of 
r; on under self-gravity can naturally enhance the ambi- 
ent density to high values. The self-g ravity of gas was ig- 
nored in the numerical simulations of Milosavl jevic et al. 
( |2009b| and jPark fc Ricotti| ( |20lT| ), on the ground that 
the BH's gravity dominate at r < r- lon . However, even 
if we start at a relatively low density initially, as long 
as the radius of the gas sphere is much larger than r lon , 
the densities just outside r- lon may reach very high lev- 
els shortly after one free-fall time scale. In this work 
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we explore this possibility and the resulting accretion 
rates with one-dimensional hydrodynamic simulations in 
spherical geometry. 

2.2. Numerical Methods 

As discussed above, there are two major feedback pro- 
cesses that affect the BH accretion, namely the photo- 
ionization heating and the radiation pressure. To study 
these processes and their effects on the BH accretion, one 
must resolve the spatial scale smaller than the sonic ra- 
dius, which is of the order of 10 14 cm for a 100 M BH. 
Moreover, in order to investigate the dependence of BH 
accretion on gas properties, a variety of gas density must 
be considered. 

In order to achieve a sufficiently high spatial resolu- 
tion, and simulate a wide range of gas density efficiently, 
we use a modified versi on of the one-dimensional hydro- 
dynamics code VH-1 by Blondin & Lufkin (1993), which 
is publicly available] 
parabolic method 



VH-1 is based on the piece- wise 



olclla & Woodward 1984), and uses 



a Lagrange-Remap approach to solve the Euler equa- 
tions: 



d t (pu) 
dt(pE) 



dtP + V 

V • {pun) - 

V • (pEu -| 



(pu)-- 
- VP = 
Pu)-- 





pa 
pu- 



+ H-C, (3) 



where the primary variables are the mass density /?, gas 
pressure P, and fluid velocity u. In our implementation, 
the acceleration field a includes both gravity and radia- 
tive pressure forces, and H and C are heating and cooling 
source terms arising from radiative feedback. 

The total energy per unit mass E is the sum of the 
kinetic energy and the internal energy e: 

E= l -v? + e, (4) 

and the internal energy is related to the pressure 
through the equation of state for an ideal gas for a given 
adiabatic index 7. 



P = peh-1). 



(5) 



The modifications we made to VH-1 include two main 
aspects. The first one is the straightforward change of a 
uniform radial grid in spherical geometry to a logarithmic 
one. This is to ensure sufficiently high resolution near the 
central BH while covering a large spherical region with a 
reasonable number of grid cells. The second one involves 
a detailed treatment of radiative feedback. Specifically, 
we added the heating and cooling terms in the energy 
conservation equation, and the radiative pressure forces 
in the momentum and energy conservation equations. 

The radiative feedback is modeled by formulating the 
cooling, heating, and radiation pressure with two param- 
eters: the ionization parameter £ and gas temperature T. 
We assume that the photo-ionized region around the BH 
can be fully described by £ and T. The ionization pa- 
rameter £ is defined as: 



L(r) 



n(r)r 2 ' 

1 http:/ /wonka. physics. ncsu.edu/pub/VH-l/ 



(6) 



where L{r) is the local luminosity, and n(r) is the hydro- 
gen number density at radius r. The BH luminosity L(0) 
is assumed to have a power law spectrum, L y ~ z^ -3 / 2 , 
and the attenuation of the luminosity is described by 



dL(r) 
dr 



-4?rr pca r abs , 



(7) 



where a abs is the radiative acceleration due to photo- 
absorption, and c is the speed of light. In solving Equa- 
tion [7J we make another simplifying assumption that the 
spectral shape of L(r) is the same as the original L(0), 
so that the heating, cooling, and radiation acceleration 
at any radius can be obtained from pre-computed tables 
of the £-T grid, 



H = n 2 H(^T) 




C = n 2 C(£ 7 T) 




a r abs= n a r abs(^ T ) 






(8) 



where the scaling on the local hydrogen density is explic- 
itly factored out, and a r line is the radiation acceleration 
mainly due to resonant Lya line scattering. 

We compute these two dimensional tables u sing the 
photo-ioniz ation code CLOUDY version 08.00 (Ferland 
et al.|1998| ), assuming a primordial elemental abundance, 
and a power-law input spectrum with index —1.5 in the 
energy range of 13.6 eV to 100 keV. In the hydrodynam- 
ics update of the momentum and energy conservation 
equations, the heating, cooling, and radiation forces are 
then interpolated from these tables, in conjunction with 
the solution of Equation [7J 

In addition to photo-absorption and resonant line scat- 
tering, Thomson scattering also contributes to the ra- 
diation pressure force, which is treated separately as 
a eont = -a 9 (r)L(0)/ L edd (r), where a 9 (r) and L edd are 
the gravitational acceleration and Eddington luminosity 
of the BH, respectively. We note that the unattenuated 
BH luminosity, L(0), is used in calculating a r contl as the 
absorbed luminosity is assumed to be re-emitted at en- 
ergies below the ionization threshold, which also con- 
tributes to Thomson scattering. 

Our feedback meth od is an improvement to those 



et al. 2005 Ciotti k 



previously used (e.g., ICiotti fc Ostriker 2001 Sazonov 
Ostrikerl | 20U7[) in that iT com- 



putes the heating, cooling, and radiation pressure forces 
more accurately, and applies to a temperature range of 
1 < T < 10 s K, much larger than the previous range 

10 4 < T < 3 x 10 7 K. As demonstrated in the main 
page, our simulations prod uce similar results as thos e by 
Milosavljevic et al.| ( |2009"b| and |Park fc Ricotti| ( [2011] ) us- 
ing more sophisticated treatments on the same problems. 
Moreover, this method allows us to simulate high-density 
clouds accurately and efficiently, which is critical in our 
investigation. 

Overall, these modifications allow us to achieve an 
unprecedentedly high spatial resolution of 10 11 cm 
(10~ 7 pc), and model a large range of gas density 

10 5 — 10 11 cm~ 3 . In particular, the ability to reach densi- 
ties several orders of magnit ude higher than the limit of 
10 7 cm" 3 in previous work by Milosavljevic et ah] ( 2009b I 
and Park & Ricotti (2011 ) is the key to our new findings. 
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Table 1 

Important Parameters and Properties of Simulations Performed 



Runs 


n (cm- 3 )W 


T-min (PC) [2] 


Tmax (pc)' 3 ' 


N -jM 

1 " grid 


Mmin l5] 


Mavg l6] 


Tcycle (yr) [7) 


Hon (IO" 3 ] 


,\ 




X 1U 


1 


1ZUU 


n it; 
U. 10 


AQ 


zuu 


O.o 


BP! 


10 7 


5 x icr 5 


1 


1200 










C 


10 5 


10" 4 


1 


1000 


0.007 


0.031 


1500 


40 


D 


5 x 10 5 


IO" 4 


1 


1000 


0.01 


0.099 


720 


25 


E 


10 6 


io~ 4 


1 


1000 


0.02 


0.15 


538 


18 


F 


5 x 10 6 


IO" 4 


1 


1000 


0.06 


0.40 


295 


11 


G 


2.5 x 10 7 


5 x icr 5 


1 


1200 


0.3 


0.70 


118 


4.6 


11 


5 x 10 7 


5 x icr 5 


1 


1200 


0.6 


0.91 


76 


3.1 


1 


10 s 


5 x 10~ 5 


1 


1200 


0.85 


1.13 


36 


2.2 


J 


5 x 10 8 


10~ 6 


1 


1200 


1.5 


1.61 


33 


0.72 


K 


10 9 


icr 6 


1 


1200 


1.72 


1.72 


21 


0.47 


L 


10 10 


5 x 10~ 7 


1 


1200 


1.82 


1.82 


12 


0.17 


Ml 10 ! 


10 11 


io- 7 


icr 2 


2000 


1.94 


1.94 


2.0 


0.057 



M Ambient density of the gas sphere. 

I 2 1 Inner radius of the simulation volume. 

t 3 l Outer radius of the simulation volume. 

M Number of points for the logarithmic radial grid. 

I 5 1 Asymptotic minimum accretion rate in unit of 10 — 6 Mq yr — 1 . 

( 6 1 Asymptotic average accretion rate in unit of 10 — 6 Mq yr — 1 . 

I 7 1 Asymptotic period of the oscillation cycle. 

I s ! Ionization radius where the ionization fraction is 0.5. 

M Only run B includes self-gravity. All other runs are without self-gravity in order to exam the main dependence of the accretion process 
on the ambient density. 
[!°] At no = 10 11 cm" 3 , 
computationally feasible. 



2.3. Initial Conditions and Parameters of Simulations 

In this study, we model the radial accretion onto a 
central BH embedded in a gas sphere. We assume a BH 
mass Mbh = 100 Mq, a radiative efficiency of 0.1, and 
the gas clouds are assumed to be metal-free and uniform 
with an ambient density hq, and an initial temperature 
T = 10 4 K. Outflow and inflow boundary conditions are 
used in the inner and outer boundaries, respectively. The 
inner and outer boundaries of the sphere are chosen in 
such a way that both sonic radius and ionization radius 
are located well within the simulation region in order to 
sufficiently resolve the accretion flow. 

A total of 13 simulations with different initial hydrogen 
densities were performed in the present work. The most 
important parameters and properties of these models are 
listed in Table 1. In this table, column 1 lists the name of 
the simulations. Columns 2-5 give the numerical param- 
eters: the ambient density of the gas cloud, the inner and 
outer radius of the simulated sphere, and number of grid 
points, respectively. Columns 6-9 give the properties of 
the accretion process: the minimum and average accre- 
tion rates, the oscillation period of the accretion process, 
and the ionization radius where the ionization fraction is 
0.5, respectively. 

Run A used the same initial cond i tions as in ot her pre- 
vious work by Milosavljevic et al. (2009b) and Park & 
Ricotti (20111. It produced similar results of many as- 
pects of the accretion process, including the general in- 
termittent pattern, accretion rates, and oscillation pe- 
riod. This confirms the previous findings that radiative 
feedback strongly suppresses BH accretion. In Run B, 
self-gravity of the gas was included. We find that self- 
gravity greatly modifies the accretion flow. It helps to 
maintain a large inflow rate and increase density buildup 



outside of the ionization radius rj on , which leads to the 
shrink of ri on . When the density reaches a critical thresh- 
old, the radiative feedback becomes less effective, so the 
accretion rate can reach Eddington limit. Runs C - M 
covered a wide range of gas density of 10 5_11 cm- 3 in 
order to explore the dependence of accretion on gas den- 
sity. In these simulations, the self-gravity is deliberately 
turned off to avoid the density enhancement due to the 
collapse outside of n on . Such controlled simulations are 
not only more computationally tractable than simply al- 
lowing the gas to collapse infinitely under self-gravity, 
they also better illustrate the main dependence of the 
accretion process on the ambient density. 

3. RESULTS 
3.1. Enhanced Accretion by Self-gravity of Gas 

The BH accretion rate is dictated by the interplay be- 
tween inflow driven by gravity and external pressure and 
outflow driven by radiation pressure. Without gas self- 
gravity, the accretion exhibits a pattern of periodic os- 
cillations, as demonstrated in Figure [I] (top panel) for 
an ambient density of 10 7 cm~ 3 . The intermittency is 
due to alternating expulsion and fallback of the gas flow 
with an average period of approximately 200 years. The 
maximum rates occasionally exceed the Eddington limit 
(A^Edd ~ 2 x 10~ 6 M Q yr -1 ), but the mean remains 
roughly constant at about 25% of the Eddington rate 
over a very long timescale. These results are in good 
agreem ent with previous ones from 2-dimensional simu - 
lations flMilosavljevic et al.|2009b||Park fc Ricotti|20lT ). 
Similar intermittent behavior were also seen in other sim- 



ulations of Bond i-like (e.g., Krumholz et ah 2005; Ciotti 
fc Ostrikerj2007 ) , or rotating accretion (e.g., |Proga et al. 
2008) on different scales. 
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Figure 1. Time evolution of accretion rates onto a 100 Mq black 
hole embedded in a uniform, primordial gas sphere with ambient 
density of 10 7 cm - 3 , without (top panel) and with (bottom panel) 
self-gravity of the gas in the simulation. The dotted line in each 
panel indicates the Eddington rate MEdd ~ 2x 10 -6 Mq yr — The 
accretion behavior exhibits a periodical oscillation caused by the 
alternating inflow driven by gravity and external pressure and out- 
flow driven by radiative feedback from the BH accretion. Without 
gas self-gravity, the average accretion is about 25% of the Edding- 
ton rate, but with self-gravity, it reaches ~ 86% of the Eddington 
rate after one free-fall timescale of the gas sphere. 



However, when self-gravity of the gas is included, the 
BH accretion pattern changes dramatically, as shown in 
Figure [l] ( bottom panel) . It is clear that the minimum 
accretion rate within one oscillation cycle gradually in- 
creases. In about one free-fall timescale of the gas sphere 
(~ 1.35 x 10 4 years), the mean accretion rate reaches a 
constant at ~ 1.72 x 1CT 6 M©yr _1 , or about 86% of the 
Eddington rate. 

The reason for such a significant enhancement is that 
the BH accretion rate critically depends on the relation 
between ionization radius r- lon and sonic radius, or more 
accurately, the standard Bondi accretion radius ignor- 
ing radiative feedback rg (~ 5 x 10 15 cm for a 100 M© 
BH in a 10 4 K cloud): the accretion is strongly sup- 
pressed by radiative feedback if r ion 3> rs, but is less 
affected if r- lon < ra- The properties of the accretion 
flow change when the r- lon = r-g boundary is crossed. At 
radius smaller than rj on , radiative heating and photo- 
ionization dominate, while the BH's gravitational field 
starts to substantially modify the gas flow at radius close 
to r-Q. When r ion < the radiation feedback is not 
important at radius larger than tb except for Thomson 
scattering. The gas flow can therefore achieve a large 
inward velocity of the order of local sound speed at rs- 
On the other hand, when ri on > rs, radiation feedback 
may reduce the local accretion rate significantly, creat- 
ing a strong shocked region between ri on and r^. It is 
this region with enhanced density that helps to maintain 
a relatively large inflow rate only limited by Thomson 
scattering. Self-gravity gradually overcomes the radia- 
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Figure 2. Time evolution of BH accretion rate for different am- 
bient gas density, n = 10 7 , 10 s , 10 9 , 10 10 , and 10 11 cm" 3 , re- 
spectively. Note the accretion curves of the last four densities are 
shifted vertically for easy comparison, and the red line represents 
the mean accretion rate, M aV g, of each simulation with the actual 
value (in unit of 10 — 6 Moyr - 1 ) given next to it. The horizontal 
axis is the scaled time, where the scale factor for different density 
is given in vertical orientation at the end of the corresponding ac- 
cretion curve. In these simulations, self-gravity of the gas is delib- 
erately turned off to avoid density enhancement due to the collapse 
outside of r; on . Such controlled simulations not only better unveil 
the important dependence of the accretion process on the ambient 
density, they are also more computationally tractable than simply 
allowing the gas to collapse infinitely under self-gravity. 



tion force and increases density buildup outside of r- 10TL , 
which leads to the shrinking of r- lon to within tb . In the 
regime where r- lon < r&, the accretion rate can approach 
to Eddington limit when the density at tb reaches a crit- 
ical value. 

3.2. Dependence of Accretion on Ambient Gas Density 

An important question is then what the critical den- 
sity would be. Figure [2] shows the time evolution of ac- 
cretion rate at several aifferent densities. As the density 
increases, the accretion rate becomes constant and ap- 
proaches to the Eddington limit. 

The oscillation period also decreases as density in- 
creases, and scales as rip ' 46 , as illustrated in Figure [5] 
(top panel). When the density reaches ~ 10 9 cm" 3 , the 
oscillation is significantly damped after an initial time pe- 
riod, and the ionization radius r- lon drops to below r&, as 
shown in Figure[3] ( bottom panel) . The radiation feedback 
becomes less effective. As a result, the mean asymptotic 
accretion rate rises to about 90% of the Eddington limit. 
A density of 10 9 cm" 3 or higher may be common in col- 
lapsing clouds in massive dark matter halos in the early 
universe, as reported in high-resolution sim ulations of the 
first protostars (e.g., Yoshida et al. 2008). This density 
is much higher than those considere d in previous work 
of similar problem (e.g., Alvarez et al. 2009 Milosavl- 
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Figure 3. Dependence of the oscillation period of the BH accre- 
tion, T cyc i c (top panel), and ionization radius r lon (bottom panel), 
on the ambient density, no, of the gas sphere. The filled circles are 
simulation data, while the solid lines are fitting curves. Both T cyc i e 
and ri on decrease with density: r cyc i c = 3 X f0 5 (rao/cm -3 ) - 0,46 
years, and r ion = 1 8.3(no/cm -3 ) -0 ' 5 pc. At no > 10 <J cm~ 3 , the 
ionization radius becomes ri on < re (rs ~ 1.6 X 10 -3 pc), the radi- 
ation feedback is weakened, and the oscillation is strongly damped. 



jevic et al. |2009b| |Park fc Ricotti||2011| ). It is therefore 
not surprising that they found severely reduced accretion 
rates. 

3.3. A Universal Correlation Between Black Hole 
Accretion Rate and Ambient Gas Density 

Figure [4] shows the relationship between the BH accre- 
tion rate and ambient density, which can be described 



as M B H = M E dd 



»0/"c 







where n c is the character- 



l+n /n c ^ 

istic density at which the asymptotic accretion rate ap- 
proaches to the Eddington limit. As shown in the Figure, 
the fitting of the simulation data gives n c = 2 x 10 8 cm -3 , 
and (3 = 0.48. At densities above 2 x 10 s cm -3 , radiative 
feedback other than Thomson scattering becomes less 
effective, and the accretion reaches the Eddington limit, 
while below that, radiative heating and photo-ionization 
feedback dominates and strongly suppresses the accre- 
tion. This relationship would imply an accretion rate 
larger than the Bondi rate at densities below ~ 10 2 cm~ 3 . 
However, because both ionization and heating timescales 
decrease as density becomes smaller, they eventually be- 
come much less than the dynamical timescale of the gas. 
Under such conditions, the radiative feedback plays a 
negligible role, and the accretion becomes Bondi-like. 

More intriguingly, we repeat the simulations for a 
BH mass range of 10 2 — 10 9 M Q , and find that this 
correlation is universal over a wide range of BH mass 
and gas density, with a simple scaling relation between 
the BH mass and the critical density, ncrit-^BH = 
(2 x lO 8 /cm~ 3 )(lO 2 /M ). This scaling relation is easy 
to understand, as the hydrodynamic equations governing 
the accretion process without self-gravity can be shown 
to depend only on u Mbh, under the assumption that 




^ 10 1 r 



Figure 4. Correlation between the BH asymptotic accretion rate 
and the ambient gas density. The filled circles represent data 
from simulations with radiation feedback, while the solid line is 
the least-square fitting, which takes the following form: A/bh = 
" o/ " c ^ where n c -'vin8 



M- 



■■ 2 X 10 s cm" 



Bdd yi+n /n c 

plot shows that, at densities below 2 X 10 s cm" 



, and p = 0.48. This 
3 , radiative feedback 



dominates and strongly reduces the accretion to be sub-Eddington, 
while above that, radiative feedback becomes less effective, and the 
accretion reaches the Eddington rate. 



the density is sufficiently high so that the local gas tem- 
perature is close to the thermal equilibrium determined 
by the heating and cooling functions. That means, for 
a given BH mass, there exists a critical density, above 
which the accretion rate can reach the Eddington limit. 

This general relation between BH accretion and am- 
bient gas density has important implications and appli- 
cations in studies of BH growth. We note that Bondi 
accretion ha s been commonly used in simulations of BH 



growth (e.g 



varez et al 
2005 



et al. 



Li et al.|2 007; Johnson & Bromm 2007[JA1- 
liuU9| |Di Matteo et al.||2005| |Springel"^a 
Hop kins et al.||2006t |Di Matteo et al.||2008)|Sijack" 



2009; Di Mat teo et al.||201ip . This simplified pre- 



scription neglects the effects of radiation feedback and 
overestimates the accretion rate by up to two orders of 
magnitude at some densities below 10 8 cm~ 3 . Our re- 
sults and the fitting formula above can serve as a more 
realistic recipe for BH accretion, and can be implemented 
directly into numerical simulations. 

4. SUMMARY 

To summarize, we have presented a set of one- 
dimensional hydrodynamic simulations of the accretion 
of a black hole embedded in a primordial gas cloud, using 
the modified grid-based VH-1 code. We include not only 
important feedback processes from the accreting black 
hole, but also self-gravity of the gas. We achieved an 
unprecedentedly high spatial resolution of 10 11 cm, and 
covered a wide range of gas density of 10 5 — 10 11 cm~ 3 . 
These advantages allowed us to study the accretion pro- 
cess in regimes not explored by previous work, and unveil 
the following new findings: 

I. The accretion behavior exhibits a periodical oscil- 
lation caused by the alternating inflow driven by 
gravity and external pressure and outflow driven 
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by radiative feedback from the BH accretion. Self- 
gravity of the gas can boost the accretion rate by 
building up high gas density which weakens the 
feedback effects. Without gas self-gravity, the av- 
erage accretion is about 25% of the Eddington rate, 
but with self-gravity, it reaches ~ 86% of the Ed- 
dington rate after one free-fall timescale of the gas 
sphere. 

2. The accretion depends strongly on the ambient gas 
density. For a given black hole mass, there exists a 
critical density above which the accretion can reach 
Eddington limit. For example, for a 100 M BH, 
the critical density is ~ 10 9 cm -3 . 

3. There exists a universal correlation between black 
hole accretion rate and ambient gas density: 

M BH = MB dd (jf£^y, where n c = 2 x 

10 s cm -3 , and j3 = 0.48. This fitting formula may 
serve as a realistic recipe for BH accretion, and 
can be implemented directly into numerical simu- 
lations. 

In the rare, high-density peaks (5 — 6 a) of the cosmic 
filaments where the first most massive dark matter halos 
collapsed, the first stars were born and might evolve into 
the first stellar-mass BHs. The extremely deep gravita- 
tional potential in these regions retained an abundant gas 
supply, and induced vigorous interactions and mergers of 
protogalaxies. Strong gravitational torques removed an- 
gular momentum from the highly shocked and dense gas 
and transported it to fuel the central BHs, while trig- 
gered global starburst on a larger scale. Seed BHs in 
such a gas-rich and dynamical environment may grow 
rapidly and co-evally with host galaxies to become the 
first luminous quasars at the cosmic dawn. 
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